home *** CD-ROM | disk | FTP | other *** search
/ Encyclopedia of Graphics File Formats Companion / GFF_CD.ISO / software / mac / nihimage / nih_src.hqx / V1.54 Source / Background.p < prev    next >
Encoding:
Text File  |  1993-10-27  |  40.0 KB  |  841 lines
  1. unit Background;
  2. {************************************************************************}
  3. {*     by Michael Castle and Janice Keller                                                                                                    *}
  4. {*     University of Michigan Mental Health Research Institute (MHRI)                                                     *}
  5. {*     e-mail address: mike_castle@ub.cc.umich.edu                                                                                   *}
  6. {************************************************************************}
  7. {*     Rolling ball and rolling disk algorithms inspired by Stanley Sternberg's article, "Biomedical        *}
  8. {*  Image Processing", IEEE Computer, January 1983.  Discussions with Rork Kuick and Tom               *}
  9. {*  Ford-Holevinski also contributed to the development of the algorithms and improvements in their   *}
  10. {*  efficiency.                                                                                                                                                *}
  11. {************************************************************************}
  12.  
  13. interface
  14.  
  15.     uses
  16.         QuickDraw, Palettes, PrintTraps, globals, Utilities, Graphics, Camera, Filters;
  17.  
  18.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  19.  
  20.  
  21. implementation
  22.  
  23.     type
  24.         IntRow = array[0..MaxLine] of Integer;
  25.         BackSubKindType = (RollingHorizontalArc, RollingVerticalArc, RollingBothArcs, RollingBall);
  26.  
  27.     var
  28.         ArcTrimPer: integer;                                  {trim off percentage of each side of the rolling ball patch}
  29.         shrinkfactor: integer;                                 {shrink the image and ball by this factor before rolling ball}
  30.         BackSubKind: BackSubKindType;                {which kind of background subtraction are we doing}
  31.         IntPlotWidth: Boolean;
  32.         Intplotwidthval: integer;
  33.         BoundRect: rect;
  34.         xxcenter, yycenter: integer;                                        {center of rectangular mask used in MinIn2DMask}
  35.         xmaskmin, ymaskmin: integer;                                   {upper left corner of mask used in AvgIn2DMask}
  36.         backgroundptr, ballptr, smallimageptr: ptr;              {ptrs to background, rolling ball, shrunk image memory}
  37.         backgroundaddr, balladdr, smallimageaddr: longint;   {addrs of background, rolling ball shrunk image memory}
  38.         patchwidth: integer;                                                      {x or y dimension of the rolling ball patch}
  39.         leftroll, rightroll, toproll, bottomroll: integer;         {bounds of the shrunk image}
  40.         Aborting: boolean;
  41.  
  42.  
  43.  
  44.     function NewPtrToClearBuf (blockSize: Size): Ptr;
  45.     {This function will return a pointer of size specified and will}
  46.     {clear the memory to zeros . This is done to create an empty bit}
  47.     {map containing nothing but white bits . }
  48.  
  49.     {MOVE . L  ( SP ) + , D0  ; get Size variable from stack}
  50.     {_NewPtr , clear           ; make pointer }
  51.     {MOVE.L  A0 , ( SP )       ; return pointer }
  52.     {MOVE.W D0, MemErr     ; set up MemErr }
  53.     inline
  54.         $201F, $A31E, $2E88, $31C0, $0220;
  55.  
  56.  
  57.  
  58.     procedure RollBall;
  59. {*******************************************************************************}
  60. {*     RollBall 'rolls' a filtering object over a (shrunken) image in order to find the image's smooth continuous    *}
  61. {*  background.  For the purpose of explaining this algorithm, imagine that the 2D grayscale image has a third     *}
  62. {*  (height) dimension defined by the intensity value (0-255) at every point in the image.  The center of the      *}
  63. {*  filtering object, a patch from the top of a sphere having radius BallRadius, is moved along each scan line of     *}
  64. {*  the image so that the patch is tangent to the image at one or more points with every other point on the patch    *}
  65. {*  below the corresponding (x,y) point of the image.  Any point either on or below the patch during this process*}
  66. {*  is considered part of the background.  Shrinking the image before running this procedure is advised due to      *}
  67. {*  the fourth-degree complexity of the algorithm.  Care has been taken to avoid unnecessary operations (exp.      *}
  68. {*  multiplication inside loops) in this code.                                                                                                                *}
  69. {*******************************************************************************}
  70.         var
  71.             halfpatchwidth,                                {distance in x or y from patch center to any edge}
  72.             ptsbelowlastpatch,                           {number of points we may ignore because they were below last patch}
  73.             left, right, top, bottom,                   {}
  74.             xpt, ypt,                                           {current (x,y) point in the shrunken image}
  75.             xpt2, ypt2,                                      {current (x,y) point in the patch relative to upper left corner}
  76.             xval, yval,                                        {location in ball in shrunken image coordinates}
  77.             zdif,                                                  {difference in z (height) between point on ball and point on image}
  78.             zmin,                                                {smallest zdif for ball patch with center at current point}
  79.             zctr,                                                 {current height of the center of the sphere of which the patch is a part}
  80.             zadd: integer;                                    {height of a point on patch relative to the xy-plane of the shrunken image}
  81.             ballpt,                                               {index to chunk of memory storing the precomputed ball patch}
  82.             imgpt,                                               {index to chunk of memory storing the shrunken image}
  83.             backgrpt,                                          {index to chunk of memory storing the calculated background}
  84.             ybackgrpt,                                        {displacement to current background scan line}
  85.             ybackgrinc,                                      {distance in memory between two shrunken y-points in background}
  86.             smallimagewidth: longint;               {length of a scan line in shrunken image}
  87.             p1, p2: ptr;                                      {temporary pointers to background, ball, or small image}
  88.     begin
  89.         UpdateMeter(20, 'Finding Background...');
  90.         left := 1;
  91.         right := rightroll - leftroll - 1;
  92.         top := 1;
  93.         bottom := bottomroll - toproll - 1;
  94.         smallimagewidth := right - left + 3;
  95.         halfpatchwidth := patchwidth div 2;
  96.         ybackgrinc := shrinkfactor * (BoundRect.right - BoundRect.left);  {real dist btwn 2 adjacent (dy=1) shrunk pts}
  97.         zctr := 0;                                            {start z-center in the xy-plane}
  98.         for ypt := top to (bottom + patchwidth) do begin
  99.                 for xpt := left to (right + patchwidth) do {while patch is tangent to edges or within image...}
  100.                     begin                                           {xpt is far right edge of ball patch}
  101. {do we have to move the patch up or down to make it tangent to but not above image?...}
  102.                         zmin := 255;                              {highest could ever be 255}
  103.                         ballpt := balladdr;
  104.                         ypt2 := ypt - patchwidth;          {ypt2 is top edge of ball patch}
  105.                         imgpt := smallimageaddr + ypt2 * smallimagewidth + xpt - patchwidth;
  106.                         while ypt2 <= ypt do begin
  107.                                 xpt2 := xpt - patchwidth;      {xpt2 is far left edge of ball patch}
  108.                                 while xpt2 <= xpt do            {check every point on ball patch}
  109.                                     begin                                   {only examine points on }
  110.                                         if ((xpt2 >= left) and (xpt2 <= right) and (ypt2 >= top) and (ypt2 <= bottom)) then begin
  111.                                                 p1 := ptr(ballpt);
  112.                                                 p2 := ptr(imgpt);
  113.                                                 zdif := BAND(p2^, 255) - (zctr + BAND(p1^, 255));  {curve - circle points}
  114.                                                 if (zdif < zmin) then begin {keep most negative, since ball should always be below curve}
  115.                                                         zmin := zdif;
  116.                                                     end;
  117.                                             end;  {if xpt2,ypt2}
  118.                                         ballpt := ballpt + 1;          {step thru the ball patch memory}
  119.                                         xpt2 := xpt2 + 1;
  120.                                         imgpt := imgpt + 1;
  121.                                     end;  {while xpt2 }
  122.                                 ypt2 := ypt2 + 1;
  123.                                 imgpt := imgpt - patchwidth - 1 + smallimagewidth;
  124.                             end;  {while ypt2}
  125.                         if (zmin <> 0) then
  126.                             zctr := zctr + zmin;                {move ball up or down if we find a new minimum}
  127.                         if (zmin < 0) then
  128.                             ptsbelowlastpatch := halfpatchwidth    {ignore left half of ball patch when dz < 0}
  129.                         else
  130.                             ptsbelowlastpatch := 0;
  131. {now compare every point on ball with background,  and keep highest number}
  132.                         yval := ypt - patchwidth;
  133.                         ypt2 := 0;
  134.                         ballpt := balladdr;
  135.                         ybackgrpt := backgroundaddr + longint(yval - top + 1) * ybackgrinc;
  136.                         while ypt2 <= patchwidth do begin
  137.                                 xval := xpt - patchwidth + ptsbelowlastpatch;
  138.                                 xpt2 := ptsbelowlastpatch;
  139.                                 ballpt := ballpt + ptsbelowlastpatch;
  140.                                 backgrpt := ybackgrpt + (xval - left + 1) * shrinkfactor;
  141.                                 while xpt2 <= patchwidth do begin     {for all the points in the ball patch}
  142.                                         if ((xval >= left) and (xval <= right) and (yval >= top) and (yval <= bottom)) then begin
  143.                                                 p1 := ptr(ballpt);
  144.                                                 zadd := zctr + BAND(p1^, 255);
  145.                                                 p1 := ptr(backgrpt);
  146.                                                 if (zadd > BAND(p1^, 255)) then  {keep largest adjustment}
  147.                                                     p1^ := zadd;
  148.                                             end;
  149.                                         ballpt := ballpt + 1;
  150.                                         xval := xval + 1;
  151.                                         xpt2 := xpt2 + 1;
  152.                                         backgrpt := backgrpt + shrinkfactor;     {move to next point in x}
  153.                                     end;  {while xpt2}
  154.                                 yval := yval + 1;
  155.                                 ypt2 := ypt2 + 1;
  156.                                 ybackgrpt := ybackgrpt + ybackgrinc;       {move to next point in y}
  157.                             end;  {while ypt2}
  158.                     end;  {for xpt }
  159.                 if ((ypt mod 5) = 0) or not FasterBackgroundSubtraction then begin
  160.                         UpdateMeter(20 + ((ypt - top) * 70) div (bottom + patchwidth - top), 'Finding Background...');
  161.                         if CommandPeriod then begin
  162.                                 beep;
  163.                                 Aborting := true;
  164.                                 Exit(RollBall);
  165.                             end;
  166.                     end;
  167.             end;  {for ypt}
  168.     end;
  169.  
  170.  
  171.     function MinIn2DMask {(xmaskmin,ymaskmin: integer)}
  172.         : integer;
  173. {*******************************************************************************}
  174. {*     MinInMask finds the minimum pixel value in a shrinkfactor X shrinkfactor mask.                                           *}
  175. {*******************************************************************************}
  176.         var
  177.             i, j,                                           {loop indices to step through mask}
  178.             thispixel,                                  {value at current pixel in mask}
  179.             min,                                          {temporary minimum value in mask}
  180.             nextrowoffset: integer;             {distance in memory from end of mask in this row to beginning in next}
  181.             paddr: longint;                           {address of current mask pixel}
  182.             p: ptr;                                        {pointer to current pixel in mask}
  183.     begin
  184.         with info^ do begin
  185.                 min := 255;
  186.                 nextrowoffset := bytesperrow - shrinkfactor;
  187.                 paddr := ord4(PicBaseAddr) + longint(ymaskmin) * bytesperrow + xmaskmin;
  188.                 for j := 1 to shrinkfactor do begin
  189.                         for i := 1 to shrinkfactor do begin
  190.                                 p := ptr(paddr);
  191.                                 thispixel := BAND(p^, 255);
  192.                                 if (thispixel < min) then
  193.                                     min := thispixel;
  194.                                 paddr := paddr + 1;
  195.                             end;     {for i}
  196.                         paddr := paddr + nextrowoffset;
  197.                     end;     {for j}
  198.                 MinIn2DMask := min;
  199.             end; {with}
  200.     end;
  201.  
  202.  
  203.     procedure GetRollingBall;
  204. {******************************************************************************}
  205. {*     This procedure computes the location of each point on the rolling ball patch relative to the center of the     *}
  206. {*  sphere containing it.  The patch is located in the top half of this sphere.  The vertical axis of the sphere         *}
  207. {*  passes through the center of the patch.  The projection of the patch in the xy-plane below is a square.           *}
  208. {******************************************************************************}
  209.         var
  210.             rsquare,                                                                         {rolling ball radius squared}
  211.             xtrim,                                                                            {# of pixels trimmed off each end of ball to make patch}
  212.             xval, yval,                                                                     {x,y-values on patch relative to center of rolling ball}
  213.             smallballradius, diam,                                                  {radius and diameter of rolling ball}
  214.             temp,                                                                             {value must be >=0 to take square root}
  215.             halfpatchwidth: integer;                                                {distance in x or y from center of patch to any edge}
  216.             i,                                                                                    {index into rolling ball patch memory}
  217.             ballsize: Size;                                                                {size of rolling ball memory}
  218.             p: ptr;                                                                            {pointer to current point on the ball patch}
  219.     begin
  220.         smallballradius := ballradius div shrinkfactor;           {operate on small-sized image with small-sized ball}
  221.         if smallballradius < 1 then
  222.             smallballradius := 1;
  223.         rsquare := smallballradius * smallballradius;
  224.         diam := smallballradius * 2;
  225.         xtrim := (ArcTrimPer * diam) div 100;                      {only use a patch of the rolling ball}
  226.         patchwidth := diam - xtrim - xtrim;
  227.         halfpatchwidth := smallballradius - xtrim;                   {this is half the patch width}
  228.         ballsize := longint(patchwidth + 1) * longint(patchwidth + 1);
  229.         ballptr := NewPtrToClearBuf(ballsize);
  230.         p := ballptr;
  231.         for i := 0 to ballsize - 1 do begin
  232.                 xval := i mod (patchwidth + 1) - halfpatchwidth;
  233.                 yval := i div (patchwidth + 1) - halfpatchwidth;
  234.                 temp := rsquare - (xval * xval) - (yval * yval);
  235.                 if (temp >= 0) then
  236.                     p^ := round(sqrt(temp))
  237.                 else
  238.                     p^ := 0;
  239.                 p := ptr(ord4(p) + 1);
  240.             end;
  241.     end;
  242.  
  243.  
  244.     procedure InterpolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  245. {******************************************************************************}
  246. {*     This procedure uses bilinear interpolation to find the points in the full-scale background given the points *}
  247. {*  from the shrunken image background.  Since the shrunken background is found from an image composed of    *}
  248. {*  minima (over a sufficiently large mask), it is certain that no point in the full-scale interpolated                 *}
  249. {*  background has a higher pixel value than the corresponding point in the original image.                                  *}
  250. {******************************************************************************}
  251.         var
  252.             i, ii,                                                   {horizontal loop indices}
  253.             j, jj,                                                  {vertical loop indices}
  254.             hloc, vloc,                                          {position of current pixel in calculated background}
  255.             vinc,                                                   {memory offset from current calculated pos to current interpolated pos}
  256.             lastvalue, nextvalue: integer;           {calculated pixel values between which we are interpolating}
  257.             p,                                                        {pointer to current interpolated pixel value}
  258.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values between which we are interpolating}
  259.     begin
  260.         vloc := 0;
  261.         with BoundRect do begin
  262.                 for j := 1 to bottomroll - toproll - 1 do begin     {interpolate to find background interior}
  263.                         hloc := 0;
  264.                         vloc := vloc + shrinkfactor;
  265.                         for i := 1 to rightroll - leftroll do begin
  266.                                 hloc := hloc + shrinkfactor;
  267.                                 bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc);
  268.                                 bglastptr := ptr(ord4(bgnextptr) - shrinkfactor);
  269.                                 nextvalue := BAND(bgnextptr^, 255);
  270.                                 lastvalue := BAND(bglastptr^, 255);
  271.                                 for ii := 1 to shrinkfactor - 1 do begin     {interpolate horizontally}
  272.                                         p := ptr(ord4(bgnextptr) - ii);
  273.                                         p^ := lastvalue + (shrinkfactor - ii) * (nextvalue - lastvalue) div shrinkfactor;
  274.                                     end;     {for ii}
  275.                                 for ii := 0 to shrinkfactor - 1 do begin     {interpolate vertically}
  276.                                         bglastptr := ptr(backgroundaddr + (vloc - shrinkfactor) * longint(right - left) + hloc - ii);
  277.                                         bgnextptr := ptr(backgroundaddr + vloc * longint(right - left) + hloc - ii);
  278.                                         lastvalue := BAND(bglastptr^, 255);
  279.                                         nextvalue := BAND(bgnextptr^, 255);
  280.                                         vinc := 0;
  281.                                         for jj := 1 to shrinkfactor - 1 do begin
  282.                                                 vinc := vinc - right + left;
  283.                                                 p := ptr(ord4(bgnextptr) + vinc);
  284.                                                 p^ := lastvalue + (shrinkfactor - jj) * (nextvalue - lastvalue) div shrinkfactor;
  285.                                             end;     {for jj}
  286.                                     end;     {for ii}
  287.                             end;     {for i}
  288.                     end;     {for j}
  289.             end;   {with boundrect}
  290.     end;
  291.  
  292.  
  293.     procedure ExtrapolateBackground2D; {(leftroll, rightroll, toproll, bottomroll: integer; backgroundaddr: longint)}
  294. {******************************************************************************}
  295. {*     This procedure uses linear extrapolation to find pixel values on the top, left, right, and bottom edges of      *}
  296. {*  the background.  First it finds the top and bottom edge points by extrapolating from the edges of the                *}
  297. {*  calculated and interpolated background interior.  Then it uses the edge points on the new calculated,               *}
  298. {*  interpolated, and extrapolated background to find all of the left and right edge points.  If extrapolation yields *}
  299. {*  values below zero or above 255, then they are set to zero and 255 respectively.                                              *}
  300. {******************************************************************************}
  301.         var
  302.             ii, jj,                                                 {horizontal and vertical loop indices}
  303.             hloc, vloc,                                          {position of current pixel in calculated/interpolated background}
  304.             edgeslope,                                          {difference of last two consecutive pixel values on an edge}
  305.             pvalue,                                               {current extrapolated pixel value}
  306.             lastvalue, nextvalue: integer;           {calculated pixel values from which we are extrapolating}
  307.             p,                                                        {pointer to current extrapolated pixel value}
  308.             bglastptr, bgnextptr: ptr;                 {pointers to calculated pixel values from which we are extrapolating}
  309.     begin
  310.         with BoundRect do begin
  311.                 for hloc := shrinkfactor to shrinkfactor * (rightroll - leftroll) - 1 do begin     {extrapolate on top and bottom}
  312.                         bglastptr := ptr(backgroundaddr + shrinkfactor * longint(right - left) + hloc);
  313.                         bgnextptr := ptr(backgroundaddr + (shrinkfactor + 1) * longint(right - left) + hloc);
  314.                         lastvalue := BAND(bglastptr^, 255);
  315.                         nextvalue := BAND(bgnextptr^, 255);
  316.                         edgeslope := nextvalue - lastvalue;
  317.                         p := bglastptr;
  318.                         pvalue := lastvalue;
  319.                         for jj := 1 to shrinkfactor do begin
  320.                                 p := ptr(ord4(p) - right + left);
  321.                                 pvalue := pvalue - edgeslope;
  322.                                 if (pvalue < 0) then
  323.                                     p^ := 0
  324.                                 else if (pvalue > 255) then
  325.                                     p^ := 255
  326.                                 else
  327.                                     p^ := pvalue;
  328.                             end;     {for jj}
  329.                         bglastptr := ptr(backgroundaddr + (shrinkfactor * (bottomroll - toproll - 1) - 1) * longint(right - left) + hloc);
  330.                         bgnextptr := ptr(backgroundaddr + shrinkfactor * (bottomroll - toproll - 1) * longint(right - left) + hloc);
  331.                         lastvalue := BAND(bglastptr^, 255);
  332.                         nextvalue := BAND(bgnextptr^, 255);
  333.                         edgeslope := nextvalue - lastvalue;
  334.                         p := bgnextptr;
  335.                         pvalue := nextvalue;
  336.                         for jj := 1 to (bottom - top - 1) - shrinkfactor * (bottomroll - toproll - 1) do begin
  337.                                 p := ptr(ord4(p) + right - left);
  338.                                 pvalue := pvalue + edgeslope;
  339.                                 if (pvalue < 0) then
  340.                                     p^ := 0
  341.                                 else if (pvalue > 255) then
  342.                                     p^ := 255
  343.                                 else
  344.                                     p^ := pvalue;
  345.                             end;     {for jj}
  346.                     end;     {for hloc}
  347.                 for vloc := top to bottom - 1 do begin     {extrapolate on left and right}
  348.                         bglastptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor);
  349.                         bgnextptr := ptr(ord4(bglastptr) + 1);
  350.                         lastvalue := BAND(bglastptr^, 255);
  351.                         nextvalue := BAND(bgnextptr^, 255);
  352.                         edgeslope := nextvalue - lastvalue;
  353.                         p := bglastptr;
  354.                         pvalue := lastvalue;
  355.                         for ii := 1 to shrinkfactor do begin
  356.                                 p := ptr(ord4(p) - 1);
  357.                                 pvalue := pvalue - edgeslope;
  358.                                 if (pvalue < 0) then
  359.                                     p^ := 0
  360.                                 else if (pvalue > 255) then
  361.                                     p^ := 255
  362.                                 else
  363.                                     p^ := pvalue;
  364.                             end;     {for ii}
  365.                         bgnextptr := ptr(backgroundaddr + (vloc - top) * longint(right - left) + shrinkfactor * (rightroll - leftroll - 1) - 1);
  366.                         bglastptr := ptr(ord4(bgnextptr) - 1);
  367.                         lastvalue := BAND(bglastptr^, 255);
  368.                         nextvalue := BAND(bgnextptr^, 255);
  369.                         edgeslope := nextvalue - lastvalue;
  370.                         p := bgnextptr;
  371.                         pvalue := nextvalue;
  372.                         for ii := 1 to (right - left - 1) - shrinkfactor * (rightroll - leftroll - 1) + 1 do begin
  373.                                 p := ptr(ord4(p) + 1);
  374.                                 pvalue := pvalue + edgeslope;
  375.                                 if (pvalue < 0) then
  376.                                     p^ := 0
  377.                                 else if (pvalue > 255) then
  378.                                     p^ := 255
  379.                                 else
  380.                                     p^ := pvalue;
  381.                             end;     {for ii}
  382.                     end;     {for vloc}
  383.             end;   {with BoundRect}
  384.     end;
  385.  
  386.  
  387.     procedure SubtractBackground2D;
  388. {*****************************************************************************}
  389. {*     This procedure subtracts each pixel from the calculated/interpolated/extrapolated background from the  *}
  390. {*  corresponding pixel value in the original image.  The resulting image is stored in place of the original        *}
  391. {*  image.  Any pixel subtractions with results below zero are given the value zero.                                           *}
  392. {*****************************************************************************}
  393.         var
  394.             hloc, vloc,                                          {current pixel location in image and background}
  395.             pvalue: integer;                                 {difference at current pixel location}
  396.             offset,                                                 {offset in memory from beginning of original image to current scan line}
  397.             backgrpt: LongInt;                              {offset to current point in background}
  398.             p: ptr;                                                {temporary pointer to image or background points}
  399.             Databand: Linetype;                           {current scan line in image}
  400.             ControlKey: boolean;
  401.     begin
  402.         backgrpt := 0;
  403.         ControlKey := ControlKeyDown;
  404.         with Info^, BoundRect do begin
  405.                 for vloc := top to bottom - 1 do begin
  406.                         GetLine(0, vloc, pixelsperline, Databand);
  407.                         for hloc := left to right - 1 do begin
  408.                                 p := ptr(backgroundaddr + backgrpt);
  409.                                 pvalue := Databand[hloc] - BAND(p^, 255);
  410.                                 if ControlKey then
  411.                                     pvalue := BAND(p^, 255);
  412.                                 if pvalue < 0 then
  413.                                     Databand[hloc] := 0
  414.                                 else
  415.                                     Databand[hloc] := pvalue;
  416.                                 backgrpt := backgrpt + 1;
  417.                             end;     {for}
  418.                         offset := LongInt(vloc) * BytesPerRow;
  419.                         p := ptr(ord4(PicBaseAddr) + offset);
  420.                         BlockMove(@Databand, p, pixelsperline);
  421.                     end;  {for}
  422.             end;     {with}
  423.     end;
  424.  
  425.  
  426.     procedure Background2D;
  427. {******************************************************************************}
  428. {*     This procedure implements a rolling-ball algorithm for the removal of smooth continuous background       *}
  429. {*  from a two-dimensional gel image.  It rolls the ball (actually a square patch on the top of a sphere) on a       *}
  430. {*  low-resolution (by a factor of 'shrinkfactor' times) copy of the original image in order to increase speed     *}
  431. {*  with little loss in accuracy.  It uses interpolation and extrapolation to blow the shrunk image to full size.     *}
  432. {******************************************************************************}
  433.         var
  434.             tport: Grafptr;
  435.             i,                                     {loop index for shrunk image memory}
  436.             backgroundsize,              {size of the background memory}
  437.             smallimagesize: Size;     {size of the shrunk image memory}
  438.             p: ptr;                             {pointer to current pixel in shrunk image memory}
  439.             table: FateTable;             {not used}
  440.     begin
  441.         ShowWatch;
  442.         UpdateMeter(0, 'Building Rolling Ball...');
  443.         GetPort(tPort);
  444.         with Info^ do begin
  445.                 SetPort(GrafPtr(osPort));
  446.                 BoundRect := roiRect;
  447.             end;
  448.         GetRollingBall;                                                                  {precompute the rolling ball}
  449.         UpdateMeter(3, 'Building Rolling Ball...');
  450.         balladdr := ord4(ballptr);
  451.         with BoundRect do begin
  452.                 leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  453.                 rightroll := right div shrinkfactor - 1;                      {on which to roll ball}
  454.                 toproll := top div shrinkfactor;
  455.                 bottomroll := bottom div shrinkfactor - 1;
  456.                 backgroundsize := longint(right - left) * longint(bottom - top);
  457.                 backgroundptr := NewPtrToClearBuf(backgroundsize);
  458.                 Aborting := backgroundptr = nil;
  459.                 backgroundaddr := ord4(backgroundptr);
  460.                 smallimagesize := longint(rightroll - leftroll + 1) * longint(bottomroll - toproll + 1);
  461.                 smallimageptr := NewPtrToClearBuf(smallimagesize);
  462.                 Aborting := Aborting or (smallimageptr = nil);
  463.                 smallimageaddr := ord4(smallimageptr);
  464.                 if not aborting then begin
  465.                         UpdateMeter(6, 'Smoothing Image ');
  466.                         filter(unweightedAvg, 1, table);                                {smooth image before shrinking}
  467.                         UpdateMeter(10, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  468.                         for i := 0 to smallimagesize - 1 do begin                {create a lower resolution image for ball-rolling}
  469.                                 p := ptr(smallimageaddr + i);
  470.                                 xmaskmin := left + shrinkfactor * (i mod (rightroll - leftroll + 1));
  471.                                 ymaskmin := top + shrinkfactor * (i div (rightroll - leftroll + 1));
  472.                                 p^ := MinIn2DMask;                            {each point in small image is minimum of its neighborhood}
  473.                             end;
  474.                         if not aborting then begin
  475.                                 Undo;        {restore original unsmoothed image}
  476.                                 RollBall;
  477.                             end;
  478.                     end
  479.                 else
  480.                     beep;
  481.                 if not Aborting then begin
  482.                         UpdateMeter(90, 'Interpolating Background...');
  483.                         InterpolateBackground2D;                              {interpolate to find background interior}
  484.                         UpdateMeter(95, 'Extrapolating Background...');
  485.                         ExtrapolateBackground2D;                             {extrapolate on top and bottom}
  486.                         UpdateMeter(98, 'Subtracting Background...');
  487.                         SubtractBackground2D;                                  {subtract background from original image}
  488.                         UpdateMeter(100, 'Subtracting Background...');
  489.                     end;
  490.             end;   {with boundrect}
  491.         DisposPtr(backgroundptr);                           {free up background, rolling ball, shrunk image memory}
  492.         DisposPtr(ballptr);
  493.         DisposPtr(smallimageptr);
  494.         DisposeWindow(MeterWindow);
  495.         MeterWindow := nil;
  496.         SetPort(tPort);
  497.     end;
  498.  
  499.  
  500.     procedure RollArc (left, rightminusone, diam: integer; var background, ballpoints: IntRow; var Dataline: Linetype);
  501.         var
  502.             xpt, xpt2, xval, ydif, ymin, yctr, bpt, yadd: integer;
  503.     begin
  504.         for xpt := left to rightminusone do begin
  505.                 background[xpt] := -255;         {init background curve to minimum values}
  506.             end;
  507.         yctr := 0;                                   {start y-center at the x axis}
  508.         for xpt := left to (rightminusone + diam - 1) do {while semicircle is tangent to edges or within curve...}
  509.             begin                                       {xpt is far right edge of semi-circle}
  510. {do we have to move the circle?...}
  511.                 ymin := 256;                          {highest could ever be 255}
  512.                 bpt := 0;
  513.                 xpt2 := xpt - diam;                {xpt2 is far left edge of semi-circle}
  514.                 while xpt2 <= xpt do            {check every point on semicircle}
  515.                     begin
  516.                         if ((xpt2 >= left) and (xpt2 <= rightminusone)) then begin  {only examine points on curve}
  517.                                 ydif := dataline[xpt2] - (yctr + ballpoints[bpt]);                {curve minus circle points}
  518.                                 if (ydif < ymin) then begin {keep most negative, since ball should always be below curve}
  519.                                         ymin := ydif;
  520.                                     end;
  521.                             end;  {if xpt2 }
  522.                         bpt := bpt + 1;
  523.                         xpt2 := xpt2 + 1;
  524.                     end;  {while xpt2 }
  525.                 if (ymin <> 256) then{if we found a new minimum...}
  526.                     yctr := yctr + ymin;   {move circle up or down}
  527. {now compare every point on semi with background,  and keep highest number}
  528.                 xval := xpt - diam;
  529.                 xpt2 := 0;
  530.                 while xpt2 <= diam do begin  {for all the points in the semicircle}
  531.                         if ((xval >= left) and (xval <= rightminusone)) then begin
  532.                                 yadd := yctr + ballpoints[xpt2];
  533.                                 if (yadd > background[xval]) then  {keep largest adjustment}
  534.                                     background[xval] := yadd;
  535.                             end;
  536.                         xval := xval + 1;
  537.                         xpt2 := xpt2 + 1;
  538.                     end;  {while xpt2}
  539.             end;  {for xpt }
  540.     end;
  541.  
  542.  
  543.     function MinIn1DMask (var Databand: LineType; xcenter: integer): integer;
  544. {*******************************************************************************}
  545. {*     MinIn1DMask finds the minimum pixel value in a (2*shrinkfactor-1) mask about the point xcenter in the *}
  546. {*  current line.  This code must run FAST because it gets called OFTEN!                                                                   *}
  547. {*******************************************************************************}
  548.         var
  549.             i,                                                                              {index to pixels in the mask}
  550.             temp: integer;                                                          {temporary minimum value}
  551.     begin
  552.         temp := Databand[xcenter - shrinkfactor + 1];
  553.         for i := xcenter - shrinkfactor + 2 to xcenter + shrinkfactor - 1 do
  554.             if (Databand[i] < temp) then
  555.                 temp := Databand[i];
  556.         MinIn1DMask := temp;
  557.     end;
  558.  
  559.  
  560. {******************************************************************************}
  561. {*     This procedure computes the location of each point on the rolling arc relative to the center of the circle     *}
  562. {*  containing it.  The arc is located in the top half of this circle.  The vertical axis of the circle passes through  *}
  563. {*  the midpoint of the arc.  The projection of the arc on the x-axis below is a line segment.                                 *}
  564. {******************************************************************************}
  565.     procedure GetRollingArc (var arcpoints: IntRow; var arcwidth: integer);
  566.         var
  567.             xpt,                                                                                 {x-point along arc}
  568.             xval,                                                                               {x-point in arc array}
  569.             rsquare,                                                                         {shrunken arc radius squared}
  570.             xtrim,                                                                            {points to be trimmed off each end of arc}
  571.             smallballradius,                                                            {radius of shrunken arc which actually rolls}
  572.             diam: integer;                                                                 {diameter of shrunken arc's circle}
  573.     begin
  574.         smallballradius := ballradius div shrinkfactor;            { operate on small-sized image with small-sized ball}
  575.         rsquare := smallballradius * smallballradius;
  576.         for xpt := -smallballradius to smallballradius do        { find the ballpoints for arc based at  (x,y)=(0,0) }
  577.             begin
  578.                 xval := xpt + smallballradius;                                     {offset, can't have negative index}
  579.                 arcpoints[xval] := round(sqrt(abs(rsquare - (xpt * xpt))));  {Ys are positive, top half of circle}
  580.             end;
  581.         diam := smallballradius * 2;
  582.         xtrim := (ArcTrimPer * diam) div 100;                       {how many points to trim off each end}
  583.         arcwidth := diam - xtrim - xtrim;
  584.         for xpt := -smallballradius to smallballradius - xtrim - xtrim do begin
  585.                 xval := xpt + smallballradius;
  586.                 arcpoints[xval] := arcpoints[xval + xtrim];
  587.             end;
  588.         for xpt := smallballradius - xtrim - xtrim + 1 to smallballradius do begin
  589.                 xval := xpt + smallballradius;
  590.                 arcpoints[xval] := 0;
  591.             end;
  592.     end;
  593.  
  594.  
  595.     procedure ExtrapolateBackground1D (var Backline, Dataline: LineType; background: IntRow; leftroll, rightroll: integer);
  596. {******************************************************************************}
  597. {*     This procedure uses linear extrapolation to find pixel values on the left and right edges of the current        *}
  598. {*  line of the background.  It finds the edge points by extrapolating from the edges of the calculated and               *}
  599. {*  interpolated background interior.  If extrapolation yields values below zero or above 255, then they are set *}
  600. {*  to zero and 255 respectively.                                                                                                                               *}
  601. {******************************************************************************}
  602.         var
  603.             i,                                                                             {index to edges of background array}
  604.             hloc,                                                                       {}
  605.             edgeslope: integer;                                                 {}
  606.     begin
  607.         with BoundRect do begin
  608.                 edgeslope := (background[leftroll + 1] - background[leftroll + 2]) div shrinkfactor;
  609.                 for i := left to shrinkfactor * (leftroll + 1) - 1 do begin     {extrapolate on left edge}
  610.                         hloc := shrinkfactor * (leftroll + 1) - 1 + left - i;
  611.                         if (Backline[hloc + 1] + edgeslope < 0) then
  612.                             Backline[hloc] := 0
  613.                         else if (Backline[hloc + 1] + edgeslope > Dataline[hloc]) then
  614.                             Backline[hloc] := Dataline[hloc]
  615.                         else
  616.                             Backline[hloc] := Backline[hloc + 1] + edgeslope;
  617.                     end;
  618.                 edgeslope := (background[rightroll] - background[rightroll - 1]) div shrinkfactor;
  619.                 for hloc := shrinkfactor * rightroll to right - 1 do begin     {extrapolate on right edge}
  620.                         if (Backline[hloc - 1] + edgeslope < 0) then
  621.                             Backline[hloc] := 0
  622.                         else if (Backline[hloc - 1] + edgeslope > Dataline[hloc]) then
  623.                             Backline[hloc] := Dataline[hloc]
  624.                         else
  625.                             Backline[hloc] := Backline[hloc - 1] + edgeslope;
  626.                     end;
  627.             end;     {with}
  628.     end;
  629.  
  630.     procedure Background1D;
  631.         var
  632.             tport: Grafptr;
  633.             hloc, vloc, arcwidth, leftroll, rightroll, numpixels: integer;
  634.             left, right, top, bottom: integer;                      {image bounds; ROTATED if RollingVerticalArc}
  635.             i, j, maskwidth: integer;
  636.             background, arcpoints: IntRow;
  637.             offset: LongInt;
  638.             p: ptr;
  639.             Dataline, Backline, Smalldataline: Linetype;
  640.             str: str255;
  641.     begin
  642.         ShowWatch;
  643.         UpdateMeter(0, concat('Shrinking Image ', long2str(shrinkfactor), ' times...'));
  644.         GetPort(tPort);
  645.         with Info^ do begin
  646.                 SetPort(GrafPtr(osPort));
  647.                 BoundRect := roiRect;
  648.             end;
  649.         GetRollingArc(arcpoints, arcwidth);
  650.         maskwidth := shrinkfactor + shrinkfactor - 1;
  651.         case BackSubKind of
  652.             RollingHorizontalArc:  begin
  653.                     left := BoundRect.left;
  654.                     top := BoundRect.top;
  655.                     right := BoundRect.right;
  656.                     bottom := BoundRect.bottom;
  657.                     numpixels := Info^.pixelsperline;
  658.                     str := 'Rolling Disk Horizontally...';
  659.                 end;
  660.             RollingVerticalArc:  begin
  661.                     left := BoundRect.top;
  662.                     top := BoundRect.left;
  663.                     right := BoundRect.bottom;
  664.                     bottom := BoundRect.right;
  665.                     numpixels := Info^.nlines;
  666.                     str := 'Rolling Disk Vertically...';
  667.                 end;
  668.         end;     {case}
  669.         leftroll := left div shrinkfactor;                                  {left and right edges of shrunken image or roi}
  670.         rightroll := right div shrinkfactor - 1;                      {on which to roll arc}
  671.         for vloc := top to bottom - 1 do begin  {for ROI}
  672.                 case BackSubKind of
  673.                     RollingHorizontalArc: 
  674.                         GetLine(0, vloc, numpixels, Dataline);
  675.                     RollingVerticalArc: 
  676.                         GetColumn(vloc, 0, numpixels, Dataline);
  677.                 end;     {case}
  678.                 for i := leftroll + 1 to rightroll do
  679.                     smalldataline[i] := MinIn1DMask(Dataline, shrinkfactor * i - 1);
  680.                 RollArc(leftroll + 1, rightroll, arcwidth, background, arcpoints, smalldataline);  {roll arc on one line}
  681.                 for i := leftroll + 1 to rightroll do begin           {interpolate to find interior background points}
  682.                         hloc := shrinkfactor * i - 1;
  683.                         Backline[hloc] := background[i];
  684.                         for j := 1 to shrinkfactor - 1 do
  685.                             Backline[hloc - j] := background[i - 1] + (shrinkfactor - j) * (background[i] - background[i - 1]) div shrinkfactor;
  686.                     end;
  687.                 ExtrapolateBackground1D(Backline, Dataline, background, leftroll, rightroll);
  688.                 for i := left to right - 1 do begin                                {subtract background from current scan line}
  689.                         Dataline[i] := Dataline[i] - Backline[i];
  690.                         if Dataline[i] < 0 then
  691.                             Dataline[i] := 0;
  692.                     end;
  693.                 case BackSubKind of
  694.                     RollingHorizontalArc: 
  695.                         with Info^ do begin
  696.                                 offset := LongInt(vloc) * BytesPerRow;
  697.                                 p := ptr(ord4(PicBaseAddr) + offset);
  698.                                 BlockMove(@Dataline, p, numpixels);            {fast whole line write}
  699.                             end;
  700.                     RollingVerticalArc: 
  701.                         PutColumn(vloc, 0, numpixels, Dataline);         {slow whole column write}
  702.                 end;     {case}
  703.                 if ((vloc mod 8) = 0) and (vloc > 16) then begin
  704.                         UpdateMeter((LongInt(vloc - top) * 100) div (bottom - top - 1), str);
  705.                         if CommandPeriod then begin
  706.                                 beep;
  707.                                 Aborting := true;
  708.                                 leave;
  709.                             end;
  710.                     end;
  711.             end;
  712.         UpdateMeter(100, str);
  713.         DisposeWindow(MeterWindow);
  714.         MeterWindow := nil;
  715.         SetPort(tPort);
  716.     end;
  717.  
  718.     procedure SetUpGel;
  719.         var
  720.             frame: rect;
  721.             AutoSelectAll: boolean;
  722.             p: ptr;
  723.     begin
  724.         if NotinBounds or NotRectangular then
  725.             exit(SetUpGel);
  726.         StopDigitizing;
  727.         AutoSelectAll := not Info^.RoiShowing;
  728.         if AutoSelectAll then
  729.             SelectAll(false);
  730.         SetupUndoFromClip;
  731.         with info^ do begin
  732.                 frame := roiRect;
  733.                 if ((LutMode = GrayScale) or (LutMode = CustomGrayscale)) and (not IdentityFunction) then
  734.                     ApplyLookupTable;
  735.                 changes := true;
  736.             end;
  737.         case BackSubKind of
  738.             RollingHorizontalArc, RollingVerticalArc: 
  739.                 Background1D;    {--------------> call background subtract <-------------------}
  740.             RollingBall: 
  741.                 Background2D;
  742.             RollingBothArcs:  begin
  743.                     BackSubKind := RollingHorizontalArc;           {remove horizontal streaks}
  744.                     Background1D;
  745.                     BackSubKind := RollingVerticalArc;               {remove vertical streaks}
  746.                     if not aborting then
  747.                         Background1D;
  748.                     BackSubKind := RollingBothArcs;                   {leave BackSubKind as we found it}
  749.                 end;
  750.         end;     {case}
  751.         UpdatePicWindow;
  752.         SetUpRoiRect;
  753.         WhatToUndo := UndoFilter;
  754.         Info^.changes := true;
  755.         ShowWatch;
  756.         if AutoSelectAll then
  757.             KillRoi;
  758.         if Aborting then begin
  759.                 Undo;
  760.                 WhatToUndo := NothingToUndo;
  761.                 UpdatePicWindow;
  762.             end;
  763.     end;
  764.  
  765.  
  766.     procedure GetBallRadius;
  767.         var
  768.             SaveRadius: integer;
  769.             canceled: boolean;
  770.     begin
  771.         SaveRadius := BallRadius;
  772.         BallRadius := GetInt('Rolling BallRadius:', BallRadius, canceled);
  773.         if (BallRadius < 1) or (BallRadius > 319) or canceled then begin
  774.                 BallRadius := SaveRadius;
  775.                 if not canceled then
  776.                     beep;
  777.             end;
  778.     end;
  779.  
  780.  
  781.     procedure DoBackgroundMenuEvent (MenuItem: integer);
  782.         var
  783.             map_array: Ptr;
  784.     begin
  785.         MeterWindow := nil;
  786.         Aborting := false;
  787.         case MenuItem of
  788.             HorizontalItem, VerticalItem:  begin
  789.                     if FasterBackgroundSubtraction then begin
  790.                             ArcTrimPer := 20;
  791.                             shrinkfactor := 4;
  792.                         end
  793.                     else begin
  794.                             ArcTrimPer := 10;
  795.                             shrinkfactor := 2;
  796.                         end;
  797.                     if MenuItem = HorizontalItem then
  798.                         BackSubKind := RollingHorizontalArc
  799.                     else
  800.                         BackSubKind := RollingVerticalArc;
  801.                     SetUpGel;
  802.                 end;
  803.             Sub2DItem:  begin
  804.                     if FasterBackgroundSubtraction then begin
  805.                             if BallRadius > 15 then begin
  806.                                     ArcTrimPer := 20;     {trim 40% in x and y}
  807.                                     shrinkfactor := 8;
  808.                                 end
  809.                             else begin
  810.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  811.                                     shrinkfactor := 4;
  812.                                 end
  813.                         end
  814.                     else begin  {faster not checked}
  815.                             if BallRadius > 15 then begin
  816.                                     ArcTrimPer := 16;  {trim 32% in x and y}
  817.                                     shrinkfactor := 4;
  818.                                 end
  819.                             else begin
  820.                                     ArcTrimPer := 12;   {trim 24% in x and y}
  821.                                     ShrinkFactor := 2;
  822.                                 end
  823.                         end;
  824.                     BackSubKind := RollingBall;
  825.                     SetUpGel;
  826.                 end;
  827.             RemoveStreaksItem:  begin
  828.                     ArcTrimPer := 20;
  829.                     shrinkfactor := 4;
  830.                     BackSubKind := RollingBothArcs;
  831.                     SetUpGel;
  832.                 end;
  833.             FasterItem: 
  834.                 FasterBackgroundSubtraction := not FasterBackgroundSubtraction;
  835.             RadiusItem: 
  836.                 GetBallRadius;
  837.         end; {case}
  838.     end;
  839.  
  840.  
  841. end.